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Abstract. We consider a hydrodynamic-type system of balance equations which is closed 
! by the dynamic equation of state taking into account the effects of spatial nonlocality. Sym- 

Oh \ metry and local conservation laws of this system are studied. A system of ODEs being 

obtained via the group theory reduction of the initial system of PDEs is investigated. The 
reduced system is shown to possess a family of the homoclinic solutions. Depending on the 
values of the parameters, the homoclinic solutions describe the solitary waves of compression 
or rarefication. The solutions corresponding to the wave of compression are shown to be un- 
stable. More likely, the waves of rarefication are stable. Numerical simulations demonstrate 
that the solitary waves of rarefication moving toward each other maintain their shape after 
the interaction. 
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1 Introduction 

This paper deals with soliton-like traveling wave (TW) solutions to some nonlinear evolu- 
tionary PDEs. TW solutions play an important role in mathematical physics. They describe 
various transport phenomena, including solitary waves on a surface of water [TJ, P waves 
in lithosphere [2], electric impulse propagation in the model of nerve axon [3}[5], moving 
fronts [6j[7], patterns on surfaces of thin films [8j|9]. 

In recent decades a number of methods enabling to obtain physically meaningful exact 
solutions were developed [101112] . Thanks to the involvement of computer algebra, the im- 
plementation of these methods is of interest for an immense amount of researchers, who, 
however, rarely pay attention to the issues of stability thereof, although the latter is a very 
important feature of any solution, which is claimed to describe observable phenomena. 

The main goal of the present paper is to study stability properties and dynamical features 
of solitary wave solutions of a hydrodynamical model, taking into account effects of spatial 



1 



non-locality. Both the effects of temporal and spatial non-locality arise naturally when 
describing the long waves propagation in media with internal structure [131115] . [8]. Besides 
the stability, we also study symmetry properties and determine local conservation laws, 
cf. [IEJCI7], admitted by the system under consideration. To the best of our knowledge, the 
symmetry properties are not directly related to the stability of TW solutions, as well as of 
the other classes of symmetry-based solutions. However, the presence of, at least elementary 
Lie (or more general) symmetry is necessary for the system in question to possess the family 
of invariant solutions, which can be either obtained in closed form, or distinguished and 
analyzed by means of qualitative methods. 

The role of conservation laws seems to be more important, but the causal relationships 
between them and the stability of TW solutions is still unknown. While the equations of the 
KdV hierarchy that possess stable soliton solutions admit a countable number of conservation 
laws, the compacton-supporting Hyman-Rosenau K(m, n) equations [18] are belived to have 
only four conservation laws [19] (for m = n it was recently proved in [20]). In spite of the 
above, compactons demonstrate stable evolution and almost elastic interaction [2"TH2"3"] . 

The paper is organized as follows. In section 2 we introduce our basic system of PDEs, 
and investigate its symmetries and local conservation laws. In section 3 we give a geometric 
insight into the structure of phase space of the dynamical system obtained via the symmetry- 
based factorization of the initial system, and formulate the conditions, which guarantee that 
the dynamical system in question possesses a one-parameter family of homoclinic solutions, 
corresponding to the solitary wave regimes. In section 4 we begin to study the spectral 
stability of the solitary wave solutions, and state the conditions which guarantee the stability 
of the limiting stationary solutions. In section 5 we study the discrete spectrum of the 
linearized problem, using the Evans function. In the end of this section we give the results of 
the numerical investigations of solitary waves' collision, demonstrating the re-establishment 
of their shape after the interaction. Finally, in section 6 we discuss the results obtained, and 
give an outline of future research. 

2 Hydrodynamic model describing long waves propa- 
gation in media with internal structure 

We are going to analyze a family of TW solutions, supported by the hydrodynamic-type 
system which takes into account the non-local effects [2~Hl25]: 



Here u is the mass velocity, p is the pressure, p is the density, t is time, x is the mass 
(Lagrangian) coordinate related the commonly used (Eulerian) coordinate x e in the following 
fashion: 



u t + (3p u+1 p x + <r[p v+1 p xxx + 3(1 + v)p v p x p xx + v{\ + v)p v ~ 1 pl] = 0, 
Pt + P 2 u x = 0, 





(3 > and tr / are constant parameters. 
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First of all, note that the system ([I]) admits a third-order Hamiltonian structure 



1 \ „ , , Jln , fO 2 



() () ^ 6a(pp xx + /t^)Da, + ( Q J (a(pp xxx + 3p x p xx ) + /3pp x 

where D x denotes the so-called total x-derivative, see e.g. [16] for details. 
We want to rewrite (EEl) in the form 



p t J \SH/5pJ 



whence we readily find that modulo the terms from the kernel of J we have 

{aDl + (5D x )5H/5u = -u x , 5H/5p - 



z/ + 2 

This yields the nonlocal Hamiltonian H which can be written as 

z 2 f p h 



H = ~J {Y + J —2 dp ) dX ' 
where z is a nonlocal variable such that 

z x = (aDl + (5)- l ' 2 u. 

Note that this nonlocality bears some resemblance with the ones that arise for the celebrated 
Camassa-Holm equation. 

We were unable to find a second Hamiltonian structure for ([1]), and the existence of such 
a structure is highly unlikely because the study of soliton dynamics for ([1]) suggests that this 
system is not completely integrable for the generic values of parameters. 

It is readily checked that for a ^ the system (CTJ possesses the following seven conserved 
quantities: 



Hi = J udx, 
H 2 = / dx/p, 



H 3 = J usmh((x)dx, 

if 4 = J ucosh((x)dx, (2) 

H e = J(tu si n H(x ) + ^)/((p))d X , 
// 7 = /( fM cosh K x) + sinh K x)/(C,))rfx, 



where £ = (— /3/a) 1 / 2 , and no further conserved quantities with local densities for the generic 
values of parameters in p]) seem to exist. 

Recall that a local conservation law for (JTj) is a relation of the form 

D t R = D X S 

where R, S depend only on x, t, u, p and a finite number of x- derivatives of u and p, and D t 
denotes the total derivative w.r.t. the time t (cf. [TS1[T7| for details). 
The local conservation laws associated with (T5D read as follows: 



D t (Ri) = D x {Si), 2 = 1,..., 7, 

where 



Ri 


= u, 


Si 


R2 


= Vp, 


s 2 


Rs 


= u sinh(Cx), 


S3 


i?4 


= u cosh(£a;), 


s± 


R 5 


= tu + xj p, 


s 5 


Re 


= tu smalt x) H , 

(p 


Se 


R 7 


. , . sinhf^x) 
= tu cosh{ux ) H , 


s 7 



p v+2 _ ^l pxx + l)pV( ^ 



u + 2 

j 

cosh ( (x ) fip u+ 1 p x 

c — 

smh((x) (3 p u+1 p x 

c 



- smh((x)a (p u+1 p xx + + l)p v p x ) 

- cosh(Cx) ( x (p u+1 p xx + + l)p v pl) 



(3) 

Note that for a > it is convenient to introduce to = ((3/a) 1 ^ 2 instead of £ and choose a 
slightly different basis of conservation laws instead of ([3]) in order to have real (rather than 
purely imaginary) densities and fluxes, viz. 

'' -p^ - a (p^p xx + + l)pffQ 



Ri 


= u, 




Si 


R2 


= l/p, 




s 2 


R 3 


= u sin(wx), 




S3 




= u cos(cjx), 




s 4 


R5 


= tu + x/ p, 




Ss 


Rq 


= tu sin(ux) 


— cos(u>x)/(cup), 


Se 


Ri 


= tu cos(ux) 


+ sm(u>x)/(ojp), 


s 7 



v + 2' 

u, 

cos(cux) (3 p u+1 p x 

cos(ux)a (p u+1 p xx + \)p v p x ) 



sm(ux)cr (p u+1 p xx + (v+ l)p u p 2 x ) 

sm(ux)(3p u+1 p x 



to 



(4) 

It is also readily checked that for <r ^ the system ([!]) admits the following Lie point 
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symmetries (see e.g. [16] for a precise definition): 



Qi = d/du, 
Q 2 = d/dx, 

Q 3 = d/dt, 1 ' 

Q 4 = (u + 3)td/dt -{v + l)ud/du - 2pd/dp, 

all of which have a clear physical interpretation: the first three are translations w.r.t. inde- 
pendent variables and the dependent variable u while the last one is the scaling symmetry. 
Somewhat surprisingly, this scaling symmetry does not involve x. 

For the generic values of parameters these appear to be the only (be it Lie point or higher) 
symmetries admitted by the system ([1]). However, for special values of parameters additional 
symmetries may emerge. 

For instance, if we concentrate on the physically relevant case of v > 0, we readily 
find that for v = 1 there appears an additional nonlocal symmetry. Namely, consider the 
extended (also known as potential, see e.g. [17] for details) system which consists of ([1]) and 
the equations for the potential w of the conservation law given by the second equation of 
([1]), that is, 

w x = l/p, w t = u. (6) 
The said extended system for v = 1 possesses the following additional symmetry: 

Q 5 = -t 2 d/dt + (tu - w)d/du + tpd/dp + (t 2 u - tw)d/dw. 



3 Qualitative study of a system of ODEs, describing 
the set of traveling wave solutions for (d) 

We are going to analyze a set of traveling wave (TW) solutions, having the form 

u(t,x) = U(z), p(t, x) = R(z), z = x — st, (7) 

where s is the velocity of TW. Inserting the ansatz ([7J into the second equation of the system 
(JT]) we get, after one integration, the following quadrature: 

U(z) = Ci - 



R(z) 



where C\ is the constant of integration. In what follows, we assume that C\ = s/Ri, where 
< Ri — const. Such a choice immediately leads to the following asymptotic behavior: 



lim u(t, x) = 0, 

\z\ — > oo 



lim p{t, x) = R\. 



Upon inserting the ansatz ([7j into the first equation of the system (CD), and using the 
equation (JHJ), we obtain the second order ODE 







R + v + 2 



R 



v+2 



a 



+1 d 2 R 
du 2 



+ (v+l)R» 



~dR~ 


2" 


du 





E. 



(9) 
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where 

B ^ + ^2< + - < 10 > 

is an integration constant, determined from the conditions on infinity. 
Let us write the equation (JHD as the first-order dynamic system: 



dR _ y 

dz 
dY 
dz 



{aR^ 2 )- 1 {ER - [s 2 + £- 2 R u+z + o{y + l)R» +l Y 2 } } 



Our aim is to state the conditions assuring the existence of a one-parameter family of solitary 
wave solutions, represented by the solutions of the dynamic system (1TT1) bi-asymptotic to 
a saddle point. To begin with, we note that all the stationary points of the system (11 II) 
are located on the horizontal axis. They are determined by the solutions of the algebraic 
equation 

P(R) = -^R u+3 - ER + s 2 = 0. (12) 

It can be easily seen that one of the roots of equation (IT21) coincides with R x . The location 
of the second real positive root depends on relations between the parameters. If v + 3 > 1 
then there exists another stationary point (R2, 0), located in the positive half-plane. If s 
satisfies the inequality 

s 2 > 4, = PRi + \ (13) 

then R 2 > R\. In case when 

s 2 < 4, = mr 3 , (14) 

the second critical point satisfies the inequality R2 < Ri- 

If v is a natural nuber or zero (which we will assume later on), then the polynomial P{R) 
has the following representation: 

P(R) = (R- R 1 )(R- R 2 )^(R), 

where function ^{R) is positive, for positive R. This is also holds true for any v > —2, or, 
in other words, whenever the function R u+Z is concave in the positive half-space. 
Analysis of the linearization matrix 



M(Ri, 0) 



1 



1,2, j^i (15) 



for the system (TTTT) shows, that the stationary points Ai(Ri,0) is a saddle in the following 

two cases 

• when o > and the inequality (jTBl takes place; 

• when o < and the inequality ( THl) takes place. 

In both of the cases the stationary point A 2 (R2,0) is a center. Thus, the system ffTTj) has 
only such stationary points, which are characteristic to the Hamiltonian system. This cir- 
cumstance suggests that there could exist a Hamiltonian system equivalent to (TTTT) . Indeed, 
the following assertion holds: 
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Theorem 1 The dynamic system < f77]j is equivalent to the following Hamiltonian system 

% = 2R»(ER- [s 2 + ^R» +3 + o{y + 1)R^Y 2 ] ) = -§§, 
w/iere = 2 aRl 2u + 1 )-f, and 

al az 1 

H = 2s 2 - + P - ffly+Q + aF 2 i? 2 ^ +1 ) - 2E- . (17) 

v + 1 (z/ + 2) 2 i/ + 2 v ; 

For a > 0, the qualitative study of the system (TX6T) is performed in [23] . The main 
result obtained there can be stated as follows. If u > —2 and s 2 > f3 Ri +3 , then the system 
possesses the solution hi- asymptotic to the saddle point A\(Ri, 0). 

Now we are going to show the existence of the homoclinic solution in the case a < 0. 
Since the system (1161) is Hamiltonian, any solution thereof can be implicitly described by the 
equation H(R, Y) = C for some C = const. Thus, the saddle separatrices can be presented 
in the form 

Y = ± — y^R— (is) 



where 

Q = 2 s 2 — + R 2{v+2) - 2 E H u 

v u + 1 (u + 2) 2 z/ + 2 

and H\ = H(Ri, 0). It is seen from the formula ( fT8l) . that the incoming and outgoing saddle 
separatrices are symmetric with respect to the horizontal axis. So, it is sufficient to consider 
the separatrice placed to the left of the saddle point and belonging the upper half-plane. Let 
us show that the following statement holds. 



Lemma 2 The function Q grows with the growth of \R — Ri\ when \R — Ri\ is sufficiently 
small. 

Proof. Differentiating the function Q(R) twice we have: 

Q'(R) = 2 s 2 R v + R 2v+3 -2E R v+ \ 



Q"(R) = 2s 2 v R v ~ l + — '— (2u + 3)R 2 " +2 -2{u + l)E R v . 



Hence 



Q(R) = Q(Ri) + Q'(Ri) (R - Ri) + \ Q"(Ri) (R ~ Rif + 0{\R - Ri\ 3 ) = 

= R u i~ 1 (fiR^-s) (R- i?x) 2 + 0(\R- Ri\ 3 ). 

As (3 Ri +3 — s > 0, the function Q(R) is growing in some vicinity of the point Ri. 
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Now observe that 



Hi = 2s- 



Ri 



2s- 



+ 







> + 2) : 



+ 



- 2 



/3 



Ri +2 



il 



2(i/+2) 



2£ 



+ 







(i/ + 1)(i/ + 2) 



z/ + 2 
2(z/ + 2) 



z/ + 2 
< +2 



1? 



u+3 



So if the inequality 



iy + 1) 



> o 



2(i/ + 2) 

holds, then ifi is positive and Q(i?) becomes negative as R — > +0. Under such conditions, the 
left branch of the function K(-R), which is located at the positive half-plane, must intersect 
the horizontal axis at some point R*. In analogy with [21], it can be shown that Y(R) forms 
in this point the right angle with the horizontal axis, and the inequalities < i?* < R± take 
place. Thus, we have the following result. 



Theorem 3 //i/ > 0, and the inequalities 



R^ 3 < s 2 <(5R\ 



u+3 



(19) 



2(z/ + 2) 

hold, then the system (T7]j possesses a homoclinic solution, formed by the separatrices of the 
saddle point (Ri, 0). The homoclinic loop corresponds to the soliton-like solution of the 
source system that describes the wave of rarefication. 

Below we give the typical phase portraits of the system ([T]). To be specific, the plots 
were made for = 1.75 and v — 0. The remaining parameters varied from case to case. 
The phase portrait corresponding to a = 3, s — 1.6 and R± — 1 is shown in Figure [1] a. 
The phase portrait of the system contains a closed loop which corresponds to the wave of 
compression. 





Figure 1: Phase portraits of the dynamic system ([TT 
a) a — 3, s — 1.6, R\ — 1; b) a — —3, s = 1.6, R\ 
i?i = 1.282. 



obtained for v 
= 1.282; and c) 



a 



0, and (3 - 
—3, s 



1.75: 
= 0.7, 



Figure CD b plotted for a = — 3 < 0, s = 1.6 and R\ = 1.282 contains the closed loop 
directed towards the vertical axis. This loop corresponds to the solitary wave of rarefication. 
Figure CD c plotted for a = — 3 < 0, R\ = 1.282 and s = 0.6 does not contain the closed 
loop, since the inequalities f|T9|) are not satisfied. 



S 



4 Spectral stability of the stationary solutions 

In the study of the spectral stability of TW solutions, it is helpful to pass to new independent 
variables 

t — tj Z — S t • 

in which the invariant TW solutions ([7]) become stationary Since the main part of the 
analysis is made numerically we restrict our further consideration to the case v = 0. In the 
new variables the system (CQ) reads as follows: 



Uf - SU- Z + (3pPz + o\PPzzz + 3 PzPzz] = 0, 

Pi ~ s Pz + P 2 us = 0, 



(20) 



(for the sake of simplicity, the bars will be omitted henceforth). We restrict ourselves to the 
analysis of the spectral stability of the TW solution (U(z), R(z)), and use the standard 
ansatz 

u(t,x) = U(z) + eexp[\t]f(z), p(t, x) = R(z) + e exp[At] g(z) (21) 

where A is the spectral parameter, and |e| <C 1. 

Inserting the ansatz (|2"T|) into the system (|2"UI) . we obtain, up to the O (|e| 2 ), the linearized 

system 



fX - sf + Rag'" + 0gR! + 3ag"R' + g'((3R + 3aR") + gaR" 
R 2 f - sg' + g(X + 2RU') = 

where (■)' denotes the derivatives with respect to z. 

Let us present (1221) as the first-order dynamical system 

Y' = AY, 







(22) 



(23) 



where Y = (g, r), \, fY 



A 



( 1 \ 

10 

CL\ 0,2 O3 Q4 

\ a 5 a G / 



PR' + aR'" + R~ 2 s (A + 2RU') 



-, a 2 



s 2 - f3R 3 - 3R 2 aR" 



«3 



3R' X 
IP ° 4 " ~R^- 



a 5 



aR ' ~* aR3 

— ; , a 6 = — . Since R(z) and U(z) tend to their limiting values (i?i and 0, 

R 2 R 2 

respectively ) as \z\ increases, the dynamic system under study asymptotically tends to the 
system with constant coefficients: 

(24) 



Y — A^Y, 



where 



A r 



1 



sX s 2 - f3R\ 



aR\ <tR\ 
X s 



\ Rl 



R{ 





1 











A 



aR x 




It is obvious, that the linearized system (122]) can be treated as the spectral problem 

Ly = Xy, y=(f,g) T , (25) 

for the operator 

/ -sd z Rad zzz + 3aR'd zz + (f3R + 3oR")d z + 0R' + oK" \ 
~ V R 2 d z 2RU' -sd z ) ' 

Recall, that the set of all possible values of A G C for which the equation (1251) has nontrivial 
solutions is called the spectrum of the operator L. The homoclinic solution (U(z), R(z)) is 
said to be spectrally stable if any possible eigenvalue A belongs to the set C \J {0} of the 
complex plane |26j . 

Remark 1 It follows from the translation invariance of the system l[2U\) . that zero belongs 
to the spectrum of L. 

As usually [25], we distinguish the continuous spectrum a cont C C, and the discrete 
spectrum o^iscr C C. Being somewhat informal, we can treat a cont as the subset responsible 
for the stability of the stationary solutions 0), and o~di S cr for the stability of the solution 
(U(z), R{z)) itself. 

Now we are going to state the conditions which guarantee that a cont C C~. In the 
limiting case \z\ — > ±oo, the variational system turns into the linear system with constant 
coefficients: 

- a Ri g'" + +s f - (3 R 1 g' = Xf, (26) 
s g' - R\f = A g, 

We assume in addition that the eigenvectors (f(z), g(z)) belong to the space of tempered 
distributions S'{1Z) |30j. With this assumption, location of the continuum spectrum can be 
defined by means of the Fourier transformation. Applying the Fourier transformation to the 
system , yields: 

a) ( fe) \ _ ( A+*o, iHRi{ve-P) \ ( /(e) ^ (27) 

where f(£), g(£) are the fourier transforms of f(z), g(z), respectively. Equating the deter- 
minant of the matrix M(£, A) to zero, we obtain the expression for eigenvalues: 

A li2 = -i£a ±yf(ae-l3)eRl, £ G R. (28) 
Thus, the following statement holds. 

Statement 4 If a < 0, and (5 > 0, then the continuous spectrum does not intersect the 
positive half-plane C + . In case when a > 0, solution represented by the stationary point 
(Ri, 0) is unstable. 
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5 Numerical study of the discrete spectrum and the 
dynamical behavior of solitary wave solutions 

An efficient tools for studying discrete spectrum of linearized operator is provided by the 
Evans function E(\), which is the analytic function of the spectral parameter A. Zeroes of 
the Evans function correspond to the eigenvalues of the linearized operators, that belong to 
the discrete spectrum [3]. 

The Evans function is constructed by evolving the linearized system, dependent on A, 
starting from the points of initiation lying at — oo in the unstable invariant manifold, and 
from the points at +00 lying in the stable one. The solutions (mostly extrapolated numer- 
ically) are then calculated for the fixed values of z (usually for z — 0), and the value of 
the Wronskian at this point, determines E(X). If for for some Ao the Evans function nulli- 
fies, then the intersection of the stable and unstable manifolds is nontrivial, there exists the 
corresponding eigenvector, belonging to the Hilbert space and, thus, A € eraser- 

In practice one usually draws the real versus imaginary part of E(X) when the spectral 
parameters varies along the border of some bounded domain B of the right half-plane C + , 
and determines the number N = j — 1, where j is the number of rotations of a vector tangent 
to the parametric curve (a so called winding number [2TH29] ). If N = 0, then the analytic 
function E(X) does not contain zeroes in the domain in question |30j, Ch. XV. Since we 
are not able to integrate numerically over the unbounded region, it is necessary to become 
convinced that the Evans function is nonzero for large |A| belonging to the positive half- 
plane of the complex plane (some estimations of this sort are made in the appendix 1). So 
analyzing iV for sufficiently large B (usually B is a semicircle lying in C + ) and analyzing 
the behavior of E(X) for large |A|, one can get a hint concerning the location of eraser- 




Re E x 10" 5 

Figure 2: The real versus imaginary part of E(X) for s = 1.6, (3 = 1.75, a = —3, R\ = 1.282. 
The spectral parameter varies along the border of the half-circle with the radius b = 12, 
symmetric w.r.t. the horizontal axis and separated from the vertical axis by a small offset 
a = 0.03. 
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Figure 3: The real versus imaginary part of E(X) for s = 1.6, /3 = 1.75, a = —3, R\ = 1.282. 
The spectral parameter varies along the border of the half-circle with the radius b = 20, 
symmetric w.r.t. the horizontal axis and separated from the vertical axis by a small offset 
a = 0.03. 

Construction of the Evans function is performed as follows. For large z the spectral 
problem becomes close to the equation P3$ , having the constant coefficients. Solutions to 
this equation, which are easily calculated, form on +oo the stable manifold U + , spanned by k 
independent eigenvectors {vj} k =l of the matrix A^, corresponding to the eigenvalues Xj with 
the negative real part. On — oo solutions of the equation (jUj) form the unstable manifold 
U~ , spanned by m independent eigenvectors {vj}™ =1 of the matrix Aoo corresponding to 
the eigenvalues with the positive real parts. Construction of the Evans function becomes 
possible if these two sets of vectors are complementary, i.e. when m = n — k. 

So, initializing ( 1231) with k vectors {uj},- =1 from U + and with n — k independent vec- 
tors from U~ and solving the equation towards z = 0, we get two sets of vectors: V~ = 
Oi(0, A), V 2 (0, A), Vfc(0, A)}, and V + = {V k+1 (0, A), V k+2 (0, A),K(0, A)}, analytically de- 
pending on A. These sets have nontrivial intersection if and only if A belongs to the discrete 
spectrum eraser- Therefore the analysis of zeros of the Evans function defined as 

E[X] = det [^(0, A), V 2 (0, A), V k (0, A), ^ +1 (0, A), V k+2 (0, A), V n (0, A)] (29) 

reveals the location of discrete spectrum of the operator L. 

The construction of the Evans function is usually based on the numerical methods. An 
appropriate method is dependent on the dimensions of the invariant manifolds U + and U~ . 
In our case both the stable and unstable invariant manifolds of the matrix happen to be 
two-dimensional. Prolongation of in mult i- dimensional cases encounters the well-known 
obstacles [31] , which can be overcome by employment of the theory of exterior algebras. 

So, using the wedge product, we determine a k— form in the vector space A k (C n ) built from 
the basic elements of the vector space C n . In our case n = 4, and the two- forms belonging to 
the space A 2 (C 4 ) correspond to the invariant manifolds U + and U~ . In particular, the basic 
vectors have the form w\ — e\ A e 2 , w 2 = t\ A e^, W3 = t\ A e±, w<± = e 2 A e%, W5 = e 2 A e±, 
w 6 — e 3 A e 4- Mapping the dynamic system ( 123]) into this space, we get: 

U' = A {2) U, (30) 
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where 








1 











\ 






a 3 


a i 


1 










a 6 











1 

















04 







-a 5 














1 







-a 5 


a i 


-a 6 


a 2 


«3 / 



A (2) 



The set of eigenvalues of the matrix consists of all possible sums of eigenvalues of the 
matrix A. Therefore the manifolds are the solutions of the system ( !30l) . for which 

lim e-^U^ = t» ± , 



S-±00 



(2) 

where v + {v~) is an eigenvector of the matrix Aoo , corresponding to the eigenvalue 

with the smallest negative (largest positive) real part. The Evans function can be presented 

in terms of solutions U in the following fashion: 



E(X) = exp |- J Tt(A®)dz^ U + (z, A) A U~(z, A). 

In calculating expression (I3TI) . we employ the relation 

U + (z, A) A U~(z, A) = (U~, ZU + ) R , 
where (■) is the scalar product in the space R 6 , 



(31) 
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Below we present the steps of calculations of the Evans function. 

After fixing the values of the parameters corresponding to the appearance of the homo- 
clinic loop, we choose the starting points from which we calculate the homoclinic solution. 
The boundary conditions in ±oo are changed with the boundary conditions posed in finite 
(sufficiently large) points ±M. The values R(±M) and Y(±M) should be chosen with 
the great precision. They are obtained by means of solving the linearized system. The 
eigenvectors of the matrix A are chosen in the form 

(R(±M),Y(±M)) = (R 1 , 0)+eq, 



where q=(l, ~F^/ ^^^.f 1+4 ^ ) are the eigenvectors of the matrix A^. 

The parameter e is specified by smoothly sewing-up the solutions starting from the initial 
values z = —M z = M, correspondingly. For M = 25 e = 1.885 ■ 10~ 6 . 

Next we make the change of variables 



U = exp{fj,±z}U , 
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t=0 1=0.175 




Figure 4: Numerical simulation of system ([I]) performed with R\ = 1.0, v = 0, s = 1.6, 
(3 = 1.75 and o = 3 

that reduces the system (1231) to the following form: 

U' = (A® - ^±h)U. (32) 

We integrate the system (|32|) with corresponding initial conditions from z = +M to z = +0, 
and next from z = —M to z = — 0. The function E(X) is proportional to the S£/ + ^ , 
so zeroes of this function coincide with zeroes of the function 

E(\) = (lj-,j:u + ^ . 

The results of some implementations of the algorithm are presented below. In order to 
investigate the Evans function's behavior within the domain lying in the positive half-plane, 
the Nyquist diagrams (JR,E(\), Q E(\)) are used, enabling to fix the number of complex 
eigenvalues of the linearized operator L. We construct numerically the map C 3 A — > E(X) 
for A = a + be 2t7T1 , t G [0; 1], i — \J— 1, < a,b G R. For o > the solitary wave solutions 
are unstable (and this is confirmed by the numerical simulation shown in Fig. Hj), so we do 
not analyze this case. Some of the Nyquist diagrams obtained for o < and the values of the 
parameters for which the solitary wave solution exist are shown in Figures [2] - [3j Analysis 
shows that the winding numbers in all these cases nullify, so the corresponding domains do 
not contain the values of the spectral parameter belonging to (Jdiscr- 

We have also performed the numerical experiments, in which the Cauchy problem for the 
system is solved with the solitary wave solutions taken as the initial data. 

Numerical experiments show that solitary wave solutions obtained for o > are unstable, 
see Figure HI while those corresponding to o < are stable and evolve in a self-similar mode. 
Using the fact that the reduced system does not depend on the sign of velocity s, one is 
also able to choose as the Cauchy data a pair of solitary wave solutions separated with 
suffice spatial interval and moving toward each other. Numerical simulations show that the 
wave packs manifest the soliton behavior, maintaining their shape after the interaction, see 
Figure |5] 
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t=0 t=6 t=11.5 




Figure 5: Numerical simulation of collision of two solitary waves moving towards each other. 
Calculations are performed under the following values of the parameters: -Ri = 1.282, v = 0, 
s = ±1.6 (for the left and right perturbation correspondingly), j3 = 1.75 and a = — 3 

6 Conclusions 

Thus, in the paper analysis of the hydrodynamic-type system ([TJ) is performed, and the 
existence of the stable solitary wave solution is shown. Let us note that throughout the work 
we do not use the analytical description of the solitary wave solutions. Their existence is 
proved on the basis of the qualitative study of the corresponding dynamic system. Such study 
is more relevant and informative than the attempts to find out the analytical description for 
solitary waves, for they are rarely expressed in terms of elementary (or special) functions. 
Besides, applying the qualitative methods, one is able to cover the total set of invariant 
solutions belonging to the given family. 

The main goal of this work is the study of spectral stability of solitary wave solutions. 
A stable evolution of the solitary waves of compression corresponding to the case a > is 
rather impossible because the intersection of both of the sets u cont and eraser with C + is 
nonempty. 

In the case o < the spectrum of linearized problem occurs to be separated from the 
positive half-plane of the complex plane. These conclusions are made on the basis of the 
analytic investigations of the essential spectrum and numerical study of the Evans functions 
made for the selected values of the parameters and covering sufficiently large domain of the 
half-plane C + . Besides, they are backed by the estimation of the Evans function made for 
large |A| lying in the positive half-space. 

Let us stress in conclusion that the amount of the conserved quantities admitted by the 
system ([I]) does not depend on the sign of a so the number of conserved quantities seems 
not to be a crucial thing from the point of view of the stability properties. However it is 



15 



worth noting that the conserved densities -R3, R4, Rq and R7 are strongly dependent on the 
sign of the parameter a. 

Remarkable and somewhat unexpected result obtained in numerical experiments is the 
elastic dynamics of interaction of solitary wave solutions. Let us stress that a mere possibility 
to implement such experiments is rare for the TW solutions supported by non-integral evolu- 
tionary equations occur, as a rule, for specific values of parameters, including the magnitude 
and direction of the velocity A study of the interaction of wave packs becomes possible in 
our case due to the stability of a single solitary wave solution. Another attendant circum- 
stance is that the factorized system is invariant under the sign of the velocity of solitary 
wave, which can move in both directions. 

Let us mention in conclusion that, besides the rigorous results concerning the spectral 
stability of the TW solution, a list of unsolved problem includes the issue of more general 
stability and attractive features [34T436] . the study of the interaction of the waves under the 
wide range of paremeters' values, and special treatment of selected case v = 1, for which the 
system <Q possesses the generalized symmetry. 
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Appendix 1 



Following a common practice |33j, we shall calculate the Evans function for a system of 
ODEs with constant coefficients, approximating the linearized system (I2"3"j) . Let us use for 
A the trigonometric representation A = |A| e iip , assuming that |A| >> 1, and — | < <p < |. 
Employing the scaling transformation 



9_ 

R' 



fj=\X\ 



-1/2, 



X 



im-'x, f=\xr 1/2 f, 



and passing to the new independent variable = we obtain for |A| >> 1 approximate 
system 

( g\ (g\ 



M 



7] 

X 

\f J 



d 7] 
dr X 
\fj 
where R = R{0). 

The matrix M has four distinct eigenvalues 
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(33) 



V k = (R , P R e iak , p 3 e 2ia « 



where p = (Rq/ct 



,V4 



if + {k — 1) 7T 



aR 4 e 1 ^-^) 



k = L...4. 



tr 
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Thus, in the generic case two eigenvectors belong to the unstable invariant manifold, and two 
other to the stable one. So the Evans function for the asymptotic problem is proportional 
to 

1111 



D 



Rldet 



3 i(3ai-¥j) 



„2ia 2 



e ia S 
g 2 ia 3 
e i(3a 3 -ip) 



D 2 i 04 



(34) 



Performing simple but tedious calculations, we finally obtain that 

D = 2i R 5 e 3iip/2 (e ia2 -e iai ) (e ia3 - e iai ) (e ia4 - e iai ) . 

Analysis shows that under the above restrictions this number never vanishes, which indicates 
that the Evans function is nonzero for IAI >> 1. 



Appendix 2 

Let us consider the eigenvalues of the matrix corresponding to the small values of A. 
Characteristic equation det [A^ — /jlI] =0 can be presented as follows 

/i 4 + a 2 /i 2 + Qfi/u + a = 0, (35) 



(3 Rf - s 2 2s\ 

~R&~' ai ~ R\a' a ° ~ R\a 
eigenvalues for small A 7^ 0, we present the solutions of (1351) in the form of series 



where a 2 = — hi > a i — ~^~) a o — 111 order to analyze the behavior of the 



/i = yU + A/il + ... 

For zero-order approximation we get the equation 

H 2 o (i4 + a 2 ) = 0. 



This equation has the solution /1q' 2 = of multiplicity 2, and a pair of nonzero solutions 
/Xo 4 = ±v^— «2- For «2 < or (3R\ — s 2 > 4 are real and have different signs. For 
«2 > or f3R\ — s 2 < 0, they are pure imaginary. 

The asymptotic series corresponding to nonzero /ig' 4 has the form: 

s > 2s 2 + m\ x2j 



For ^ the expansion is different: 



where /ii = -==, /x 3 



fx = A/ii + X 3 fi 3 + 
li\R\a 



When a 2 < — Y s > \J f3R\ and A > 0, then the real roots are shifted to the right, while 
zero roots give rise to a pair of real roots having different signs. So the system (1231) has 
two-dimensional unstable invariant solution. When a 2 > — > s < ^/ j3R\, a pair of positive 
roots is created from zero ones, while the pair of a pure imaginary roots gain negative real 
part. Thus, in this case the system (T23]) has two-dimensional unstable invariant manifold. 
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